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Abstract 

Palaeoclimate archives contain information on climate variability, trends and mechanisms. 
Models are developed to explain observations and predict the response of the climate system 
to perturbations, in particular perturbations associated with the anthropogenic influence. 
Here, we review three classical frameworks of climate modelling: conceptual, simulator-based 
(including general circulation models and Earth system models of intermediate complexity), 
and statistical. The conceptual framework aims at a parsimonious representation of a given 
climate phenomenon; the simulator-based framework connects physical and biogeochemical 
principles with phenomena at different spatial and temporal scales; and statistical mod- 
elling is a framework for inference from observations, given hypotheses on systematic and 
random effects. Recently, solutions have been proposed in the literature to combine these 
frameworks, and new concepts have emerged: the emulator (a statistical, computing efficient 
surrogate for the simulator) and the discrepancy, which is a statistical representation of the 
difference between the simulator and the real phenomenon. These concepts are explained, 
with references to implementations for both time-slices and dynamical applications, 
highlights: 

1. Clarify three palaeoclimate modelling frameworks 

2. Different definitions and applications of the conceptual models 

3. Comment on usage and interest of numerical simulators 

4. Show how statistical methods apply to process-based modelling 

5. Review the current use of emulators and discrepancy in palaeoclimate modelling. 

Keywords: modelling, conceptual, general circulation models, time-space process, 
inference, Bayesian palaeoclimate reconstructions 



1. Introduction 

Palaeoclimatology, as any natural science, depends on observations to generate scientific 
activity. Questions about ice ages, Dansgaard-Oeschger events, fluctuations in CO2 atmo- 
spheric concentration, would not have been raised if these phenomena had not been observed 
in the first place. Scientists attempt to model these phenomena in order to explain them, to 
predict them, raise new questions and suggest additional experiments. 
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In climate science modellers are confronted to a dilemma about model complexity. Al- 
though this dilemma is often presented in terms of a compromise about computing power, 
there is a more fundamental dichotomy. 

On the one hand, the explanatory power of a model tends to be better accepted if 
the number of ad hoc hypotheses is as small as possible, even if this comes at the cost 
of not reproducing all the details of the phenomenon. This is the celebrated Ockham's 
razor principle, that encourages us to explain as much as possible with as little as possible. 
Consistent with this line of reasoning, Held (2005) recommends systematic research on simple 
hydrodynamical models such as the dry baroclinic atmosphere. As shall be reviewed here, 
even simpler systems, based on a small number of equations, have non-trivial properties which 
are helpful to interpret glacial-interglacial cycles and abrupt events. Furthermore, small 
mathematical models are easy to communicate so that the results can be well reproduced 
by peer scientists on the basis of information available in the publication. 

On the other hand, climate scientists have developed large numerical models, which en- 
capsulate available knowledge on a huge variety of climate processes at different spatial and 
temporal scales, ranging from cloud formation to sediment kinematics. Here we will term call 
these models: simulators. Current climate simulators are impressively successful at repro- 
ducing many features of our climate system , and the idea is to use them as experimental 
substitutes to real climate system. However, running simulators is a technically involved 
operation, and it is difficult if not impossible to fully appreciate the consequences of all the 
technical choices and physical hypotheses that they contain (Winsberg, 1999). 

A number of concepts and theories may help us to articulate these different modelling 
frameworks. The starting point is that mathematical models generate uncertain informa- 
tion about the climate system, and the purpose of any climate modeler is to connect this 
information with the real world. The branch of mathematics concerned with inference in 
presence of uncertain information is nothing but statistics. 

The purpose of the present article is to review the traditional approaches to palaeoclimate 
modelling, and, then, to show how statistics may help us to progress in the different problems 
involved in the prediction and explanation of (palaeo-) climate phenomena. The diagram on 
Figure 1 may be used as a road map. The different colors constitute different frameworks of 
palaeoclimate modelling: conceptual models, simulators (which includes general circulation 
models) and statistical models. They are connected to a number of concepts, which are 
detailed in Section 2. The light gray nodes represent traditional connections between the 
modelling frameworks, while the dark grey ones refer to more recent concepts, introduced 
and commented on in Section 3. 



a list of references is available at http://www-pcmdi.llnl.gov/ipcc/subproject_publications.php. 
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2. Traditional frameworks of palaeoclimate modelling 

2.1. Conceptual models 

2.1.1. Definition and construction strategy 

In general, a 'conceptual model' is a drastically simplified representation of a complex 
process. However, even in this context, the word 'model' is used diversely in the literature 2 . 

It may be a schematic description of the mechanisms involved in a complex process, 
often supported by a visual diagram. The 'SPECMAP' model (Imbrie et ai, 1992, 1993) 
is a description of the chain of responses of the different components of the Earth climate 
system involved in glacial-interglacial cycles (Figure 2). Likewise, the results of sensitivity 
experiments with general circulation models are often summarised in the form of a diagram. 
For example, Zhao et al. (2005) summarise a mechanistic interpretation of the effect of the 
ocean on the African and Indian monsoons with two diagrams outlining the roles of trade 
winds, evaporation and stratification. 

The word 'model' may also refer to a perceptual analogy. For example, the thermohaline 
circulation was depicted as a 'conveyor belt' (Broecker, 1991) 3 , a 'leaky funnel' (Mouchet 
and Deleersnijder, 2008) and, before this, it was compared to a device that can be set 
up in the laboratory with two vessels and a system of stirrers, porous walls and capillary 
pipes materialising the effects of mixing, diffusion and advection (Stommel, 1961) (Figure 3). 
The different perceptual analogies convey different impressions about the roles of advection 
and diffusion in the ocean circulation physics. Electronic circuits are another basis for 
perceptual analogies. Hansen et al. (1984) introduced the concept of "gain", originally 
used in electronics, to quantify the effects of "feedbacks" of climate system components 
to a radiative perturbation. This electronic analogy was made more explicit by Schwartz 
(2011), and it is further discussed later in this section (Figure 4). More abstract perceptual 
analogies have been suggested, in particular the relaxation model, in which it is assumed 
that the system is being relaxed (like a spring) to various states determined by the forcing 
and the history of the system. This analogy has been much used by Paillard (Paillard and 
Labeyrie, 1994; Paillard, 1998, 2001; Paillard and Parrenin, 2004). 

Finally, 'conceptual model' may refer to a simple system of mathematical equations. 
There are different strategies exist to determine these equations: 



2 cf. also Winsberg (1999) for more discussion on the meaning given to the word 'model' in the context of 
simulation 

3 The conveyor belt logo was introduced in the November issue of the Natural History magazine and the 
(Broecker , 1991) reference quoted here is a discussion of the relevance and limit of this analogy. 



Figure 1 (following page): Conceptual representation of the three main frameworks of palaeoclimate mod- 
elling and bridging concepts. The main frameworks are represented by the three main nodes: simulators 
in red, conceptual models in blue and statistical models in green. Light grey nodes represent traditional 
connections between the main modelling frameworks, while the dark gray ones indicate concepts more re- 
cently introduced in the palaeoclimate modelling literature. Numbers in nodes refer to section numbers of 
the present article. 
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Figure 2: The Imbrie et al. (1992) conceptual model of the linear climate system response to the astronomical 
forcings associated to precession and obliquity. The model describes a sequence of initial, early and late 
responses noted S\ to $4, to the astronomical forcing (Q), each response lagging the previous one by 0.5 to 
4 ky (1 ky=l,000 years). Reproduced with permission from American Geophysical Union. 

• They may be written to translate a perceptual analogy in order to analyse and discuss 
its consequences on the system dynamics. For example, the Stommel (1961) model 
can be expressed as two equations constrained by the laws of conservation of tracers 
combined with a parameterisation of the response of the inter- vessel flow to the density 
difference between the two vessels. One non-trivial consequence of this model is that 
it features two stable states, which have been interpreted as the on and off modes 
of the thermohaline circulation. Possible implications for the future of our climate 
are discussed in Rahmstorf (2000). Stommel's model is the simplest example of a 
wider class of mathematical models called 'box models'. The Winton and Sarachik 
(1993) model, for example, was used to study ocean internal oscillations involved in 
Dansgaard-Oeschger events (Schulz et al, 2002). The "Multi-box model" (MBM) 
(Munhoven, 2007) includes 10 boxes for the ocean and was used to study climate- 
ocean-sediment interactions over the latest glacial-interglacial cycle. The BICYCLE 
model, based in part on MBM, was used to interpret isotopic records (Kohler et al, 
2005). It includes 6 ocean boxes plus a number of reservoirs for the terrestrial biosphere, 
rocks, sediments and atmosphere. 

• Another strategy for developing mathematical conceptual models consists in starting 
from fluid dynamics equations and simplify them mathematically as much as possible, 
using a procedure called 'truncation'. The technique is well-established (Saltzman, 
1962) and it was applied to continental ice flow dynamics to study ice ages (Oerlemans, 
1980; Ghil and Le Treut, 1981) and Heinrich events (Paillard, 1995). 

• Finally, the development of a conceptual model may be more heuristic, combining 
physical arguments, information obtained from sensitivity experiments with general 
circulation models and hypotheses on non-linear effects. This approach was champi- 



oned by Saltzman (2001). 



Mathematical models may then be distinguished according to their mathematical proper- 
ties, in particular, whether they are linear or non-linear, and whether they are deterministic 
or stochastic. This aspect is now further developed. 
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Figure 3: Mechanical representation of the classical three-box Stommel (1961) model, with an explicit 
representation of diffusive phenomena (through porous walls) and advective flows (through the capillary 
pipe). The Stommel model is a prototype of more complex models used nowadays (see text). Reproduced 
from the original publication. 

2.1.2. Linear models 

Consider an electronic circuit with resistors and a capacitor (Figure 4). Resistors dissipate 
current proportionally to the difference in voltage to which it is connected (the coefficient is 
the resistance)] and the voltage across the capacitor increases proportionally to the charge 
being stored on it (the coefficient is the capacity). The system is said linear if the resistances 
and capacitance are constant. 

In an analogy with the climate system, the current may viewed as the input radiation 
(in a global warming experiment) or the net accumulation of snow on the ice sheets (in 
a Milankovitch forcing experiment). The capacitor accumulates charge over time. At the 
Milankovitch scale, ice sheets play this role, since they accumulate snow mass imbalance 
over several thousands of years. The growth of ice sheets causes a tension (voltage) on the 
system, and the role of resistors is to dissipate this tension with a discharge current, so that 
tension does not grow to infinity. At Milankovitch time scales, the discharge of ice towards 
the ablation region and the oceans plays this role. In electronics it is also possible to design 
'negative resistance circuits' (Linvill, 1953). Technically, this involves operational amplifiers. 
They inject current proportionally to the voltage. Consequently, they amplify the effects of 
the forcing. In palaeoclimates, the ice albedo feedback or the boreal vegetation feedback 
may be modelled as negative resistances because they amplify the forcing. The combination 
of positive and negative resistances may be summarised by an equivalent 'net resistance', as 
indicated on Figure 4. 

The capacitor is a crucial component of the system. It accumulates tension, like ice sheets 
accumulate mass. The effect of the capacitor on the system dynamics is to introduce a phase 
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Figure 4: Linear electronic circuit as a metaphor for the climate system. The forcing is the input current 
Q. Depending on the framework, it may represent the annual mean radiative forcing (for a global warming 
experiment) or the net mass balance over ice sheets (for a Milankovitch forcing experiment). The capacitor 
represents the component of the climate system that effectively integrates the forcings (ocean heat content, 
or ice sheets mass) and the different resistances represent what are generally referred to as "feedbacks" in the 
climate system. The values of theses resistances are consistent with the notation introduced by Hansen et al. 
(1984): Positive resistance (A > in Hansen et al. (1984)) correspond to "negative feedbacks" (they work 
against the forcing), negative resistances, in practice obtained with negative impedance circuits, are "positive 
feedbacks" . The system output (for example: temperature) is the voltage. At a given forcing frequency, the 
output lags the input with a phase tan -1 (wCA), where ui is the angular velocity of the forcing. 

lag between the forcing (input current) and the response (output voltage), so that the forcing 
and the response may be unambiguously identified. The forcing-response phase lag is used 
by Gregory et al. (2004) to estimate the net 'feedback' factor (equivalent to the inverse of 
the net resistance in our diagram) in global warming experiments with general circulation 
model simulations. In palaeoclimate applications, the phase lag is estimated from the data 
to distinguish the response from the forcing (Imbrie et al, 1992; Shackleton, 2000; Lisiecki 
et al, 2008). 

The assumption of linearity has three consequences. First, given that resistances are 
constant, the net resistance of the system is constant, too. Consequently, it must be positive 
otherwise tension (response) would run away to infinity. The implication is that small per- 
turbations to the input (for example: noisy fluctuations) are damped with time rather than 
being amplified. The system remains thus predictable even in presence of small fluctuations. 
Second, the output voltage contains exactly the same frequencies as the input current, and 
every frequency can be analysed separately from the others. The technique introduced in 
the SPECMAP project (Imbrie et al., 1992) is compatible with this linear framework. It 
consists in filtering climatic signals to separate them into components related to precession 
(19-24 ky band, 1 ky=l,000 years), obliquity (40-ky band) and eccentricity (100-ky band), 
and then inspect leads and lags between the astronomical forcing and the different climatic 
components, such as CO2, ocean temperatures, dust, etc. Visually, the picture proposed by 
SPECMAP is a chain linking the forcing to different components, of which the response has 
an increasing phase lag with the forcing as one goes down the response chain (Figure 2). 
This analysis approach is still used nowadays (e.g. Lisiecki et ai, 2008). 
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Third, the linear formalism cannot adequately account for the presence of saw-tooth- 
shaped 100-ky ice age cycles. At least three answers to this problem may be found in the 
literature: 

1. The first solution consists in postulating the existence of some unidentified 100-ky 
cycle that is effectively treated as a forcing in the linear framework, even if its origin is 
internal to the climate system. This is the solution proposed in the SPECMAP model 
(Imbrie et al, 1993). 

2. The second approach assumes that the effect of the forcing is non linear. For example, 
in Imbrie and Imbrie (1980), the response to insolation anomaly is greater if this 
anomaly is positive than if it is negative. The consequence of this asymmetry is that 
the input signal is at least partly rectified. In particular, the precession forcing, which 
is the product of eccentricity and the sine of the longitude of the perihelion (e x 
sin to) gives rise to a pure eccentricity signal in the output. The eccentricity signal 
contains the frequencies of 412, 94, 123 and 100 ky (Berger, 1978), and therefore these 
frequencies are also found in the output signal. Ruddiman (2006) proposed a conceptual 
model which, although not formulated mathematically, fits this category. The problem 
with this second approach is that the internal dynamics of the system remain linear. 
Only the effect of the forcing is being transformed non-linearly. Consequently, as 
we discussed above, there is no possibility of a runaway feedback phenomenon. This 
causes a problem to explain the deglaciation which occurred 400,000 years ago because 
insolation variations were weak at that time. If the system is stable, it cannot react 
dramatically to a small forcing. The problem was fully appreciated by Imbrie and 
Imbrie (1980). 

3. The third solution is to explain the emergence of an 'internal 100-ky cycle' as a result 
of non-linear interactions with the climate system. This is the subject of the next 
section. 

2.1.3. Non-linear, deterministic models 

A non-linear electronic system is one in which the effective resistance of the components 
vary with the tension on the system. This framework allows for a variety of new phenomena. 
In "some basic experiments with a vertically-integrated ice sheet model", Oerlemans (1981) 
noted one, two or even three stable solutions may co-exist, depending on the forcing. This 
configuration, typical of non-linear system dynamics, may cause transient behaviours that 
will look very different to the forcing functions. In particular, with some modest adjust- 
ments to the model equations and parameters (as in Pollard (1983)), such a model may 
exhibit self-sustained oscillations, even if the forcing is constant. In non-linear dynamics 
theory, such a system is called an oscillator. In the Pollard (1983) model, the oscillation 
arises from interactions between the lithosphere and the ice-sheet climate system, but only 
when model parameters lie in a fairly narrow range, that is suitable to represent the West 
Antarctic ice sheet but not the northern ice sheets. Consequently, investigators have been 
looking for other sources of oscillations in the climate system to explain ice ages. Saltzman 
and Maasch (1990), for example, proposed a system in which CO2 reacts non-linearly to 
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changes in ocean temperature, and ocean temperature is controlled by CO2 and continental 
ice volume. Contrary to the SPECMAP model, which is graphically represented as a linear 
chain of responses (Figure 2), the Saltzman-Maasch model is best graphically represented 
as a network with cyclic interactions (Figure 5). Some of these interactions act as stabil- 
ising factors (negative feedbacks, noted ' — '), but others may be stabilising or destabilising 
depending on the system state. This latter characteristic allows, in the Saltzman-Maasch 
model, ice ages to occur even in absence of astronomical forcing. Consequently, in a non- 
linear system, the phase relationship between two components (for example, astronomical 
forcing and ice volume) is no longer a reliable indicator of the nature of the forcing (see 
also Ganopolski and Roche (2009)). Compared to linear systems, the analysis focuses less 
on the analysis of phase lead and lag, and is more concerned with identification of stability, 
bifurcation, of synchronisation, and predictability, which are now briefly reviewed. 

Stability: The analysis of a non-linear dynamical system consists in identifying its stable 
and unstable points, that is, the states to which the system may be attracted (negative 
feedbacks dominate around this point) and the states from which it will be repelled 
(positive feedbacks dominate around this point). For example, Brovkin et al. (1998, 
2003) estimated that non-linear interaction between vegetation and the atmosphere 
may explain the co-existence of several stable states in the Sahel but not at the northern 
high latitudes; this work was inspired in part from the pioneering contribution by Ghil 
(1976), interested in the more general consequences of albedo feedbacks on the system 
stability. There may also be no stable point, which leaves us with several possibilities 
4 : (i) the system exhibits a stable self-sustained periodic oscillation; (ii) its behaviour 
is quasi-periodic (combination of several periods); (iii) it is aperiodic, with a broad 
power spectrum and complex phase-diagram figures (for example: the behaviour of 
ENSO modelled in Tziperman et al. (1994)). 

Bifurcation: A bifurcation is defined as a change in system behaviour obtained when one 
of the system parameters crosses a threshold, called the bifurcation point. Saltzman 
and Maasch (1990), for example, proposed models where the Middle Pleistocene Rev- 
olution is interpreted as a bifurcation. The parameter controlling this bifurcation is 
the background CO2 level, which is driven by tectonics. This bifurcation corresponds 
to a switch from a linear response mode to astronomical forcing (before the Middle 
Pleistocene Revolution) to a regime of non-linear synchronisation (after the Middle 
Pleistocene Revolution). This model is compatible with the findings of Lisiecki and 
Raymo (2007) who show, based on time series analysis, a transition from linear to 
non-linear regime before the Middle Pleistocene Revolution. At a different timescale, 
the transition between an 'interglacial' and a 'glacial' phase may also be interpreted 
in terms of bifurcation theory (Ditlevsen, 2009; Livina et al, 2011) 5 . Recent research 



4 see, among others, Ghil and Childress (1987); Saltzman (2001) and Crucifix (2012) for introductory 
texts on this in the palaeoclimate context 

5 The difference between the Ditlevsen (2009) and Saltzman and Maasch (1990) frameworks is quite 
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activity has focused on the search for 'early warning precursors' of a bifurcation Dakos 
et al. (2008); Scheffer et al. (2009); Ditlevsen and Johnsen (2010). 

Synchronisation: Synchronisation is defined as the phenomenon by which the natural os- 
cillation period of a system adjusts itself on the period of an external factor (Pikovski 
et al, 2001). Conceptual models such as those by Saltzman and Maasch (1990), Pail- 
lard and Parrenin (2004) and Tziperman et al. (2006) exhibit a phenomenon of syn- 
chronisation to the astronomical forcing 6 The concept of synchronisation is also useful 
to explain phenomena associated with Dansgaard-Oeschger and Heinrich events [Schulz 
et al., 2002). A review is available in Crucifix (2012). 

Predictability: Given that positive feedbacks may dominate at least at certain times in a 
non-linear systems, there is the possibility small perturbations be amplified so much 
that the exact history of the system from given initial conditions is in practice unpre- 
dictable. Some conceptual models of ice ages have this property. More specifically, 
climate perturbations occurring at strategic times may be amplified, and inflect the 
course of climate, by hastening or delaying significantly a glacial inception or a deglacia- 
tion [De Saedeleer et al, 2012). 

In summary, non-linear deterministic models widen the scope of plausible ice age theories 
and address questions at a different level than linear systems. The idea of a 'chronological 
chain of response mechanisms' is ambiguous in such systems because self-sustained loops are 
possible. Non-linear dynamical system theory provides a suitable vocabulary to study the 
effects of a network of interactions between different components of the climate system and 
the role of astronomical forcing. 

2.1.4- Stochastic models 

A stochastic process is a form of noise, that is, the succession of random numbers with 
well-defined properties such as auto-correlation and amplitude. There are two possible rea- 
sons for introducing a stochastic process in a model. One is to account for the physical 
phenomena that occur at smaller spatial or temporal scales than those that are explicitly 
resolved by the model. For example, Hasselmann (1976) introduced a stochastic process to 
account for 'weather' in a 'climate' model. The properties of the stochastic process may 
then be determined by reasoning on the physics of the unresolved process (Penland, 2007). 
For representing weather, it was proposed to deduce noise properties from quasi-geostrophic 



fundamental but its discussion is beyond the scope of the present review. 

6 The Southern Ocean foraminifera Fourier analysis by Hays et al. (1976) is often referred to as a pioneering 
demonstration of the influence of the astronomical forcing on climate. Retrospectively, it is tempting to 
interpret the title of that paper: "Variations in the Earth's orbit, pacemaker of ice ages" as a visionary 
reference to the concept of synchronisation. Hays et al. (1976) indeed perfectly realised that occurrence of 
100-ka ice ages does not fit a linear response framework of precession and obliquity. Though, as correctly 
observed by one reviewer of the present paper, their interpretation is closer in nature to the Imbrie and 
Imbrie (1980) model than of the self-sustained oscillation model. 
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(a) Saltzman and Maash (1991) model 
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Figure 5: (a) Graphical network representation of the Saltzman and Maasch (1990) model of ice ages. 
Arrows indicate causal relationships. The astronomical forcing causes variations in ice volume, which then 
are involved in a network of interactions with deep ocean temperature and CO2. Contrary to the linear 
framework, certain feedbacks may vary and change sign depending on the system state, which induces in the 
Saltzman and Maasch (1990) model a phenomenon of self-sustained oscillations. Ice ages are, in this model, 
self-sustained oscillations synchronised on the astronomical forcing. The result of a simulation reproducing 
Figure 5 of the original publication is shown in panel (b). Additive stochastic effects may randomly alter 
the timing and amplitude of ice ages. One example of a stochastic realisation is shown on panel (c). 
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theory (Majda et al, 2009). Stochastic parameterisations may also be introduced to ac- 
count for uncertainties and model errors, in order to more realistically assess the horizon 
of predictability of the model in the presence of such errors. Technically, accounting for 
stochastic effects transforms the dynamical system into a 'stochastic dynamical system', of 
which the numerical resolution is both theoretically and computationally more involved than 
that of a deterministic system (Kloeden and Platen, 1999). Stochastic processes also induce 
a number of interesting mathematical phenomena, which in turn provide us with a series of 
new concepts that may be adequate to interpret the history of climate. One such concept 
is the 'stochastic resonance', in which the stochastic process is combined with a periodic 
forcing to induce periodic oscillations in the system. Ice ages have been interpreted in the 
past as a phenomenon of stochastic resonance (Benzi et al., 1982). The concept of stochas- 
tic resonance itself has evolved and models showing different forms of stochastic resonance 
have been proposed to understand abrupt climate events (Ganopolski and Rahmstorf , 2002; 
Timmermann et al., 2003; Braun et al., 2007). 

2.2. Climate simulators 

Simulators are complex numerical systems designed to account explicitly for a large num- 
ber of different processes involved in the dynamics of the climate system. A particularity 
of simulators is that the number of variables is several order of magnitudes greater than 
the number of system parameters Claussen et al. (2002). This criteria is met by general 
circulation models as those used, for example, in the Palaeoclimate Modelling Intercompar- 
ison Project (Braconnot et al, 2007), but the definition applies also to most 'Earth System 
Models of Intermediate Complexity (EMICS)' (Claussen et al., 2002), such as, among others, 
FAMOUS (Smith, 2012), LOVECLIM (Goosse et al, 2010), CLIMBER-2 (Petoukhov et al, 
2000; Calov et al, 2005), CLIMBER-3a (Montoya et al, 2005), BERN3D (Ritz et al, 2010) 
and GENIE- 1 (Ridgwell et al, 2007). 

The word 'simulator' is uncommon in the climate literature. It is adopted here on purpose 
to distinguish it from statistical and conceptual models. Compared to the conceptual models, 
the purpose of the simulator is more general, in the sense that it is developed with many 
possible applications in mind (a same climate simulator may be used to study El-Nino and 
sea-ice variability; as another example, the GENIE-1 model (Ridgwell et al, 2007) contains 
the definition of 49 dissolved tracers, which can be enabled or disabled depending on the 
modeller's purpose). At the same time, experimental designs are quite specific: among 
others, the shapes and orography of continents, vegetation types, soil properties and ocean 
bathymetry have to be specified on a given grid and adapted for a particular climate era, 
for example the pre-industrial era or the Last Glacial Maximum. 

The interest of simulators is to generate a self-consistent picture of regional and planetary- 
scale phenomena that is compatible with physical and biogeochemical principles implemented 
at the scale of a grid-cell. A large fraction of computing resources is used to solve fluid 
dynamics equations that determine the movements of oceans, atmosphere, and possibly 
sea-ice and continental ice. The choice of spatial resolution conditions the spatio-temporal 
spectrum that can be studied with the simulator, ranging from weather forecast (hi-resolution 
models) to multi-millennial phenomena (EMICS). Physical and biogeochemical processes are 
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then introduced consistently. For example, a module of ocean sediment diagenesis finds its 
place in a simulator designed to study multi-millennial climate evolution (Brovkin et al, 
2012), but not in a weather forecast simulator. 

The simulator also involves a number of 'sub-grid parameterisations'. These are models 
representing the effects of physical or biogeochemical processes that are not explicitly re- 
solved by the simulator. The development of these parameterisations obeys the same rules 
as that of the conceptual models outlined above, and they may be diagnostic (assume an 
instantaneous equilibrium response), prognostic (the parameterisation is a dynamical sys- 
tem, for example, the prognostic cloud scheme of the UK Met Office model (Wilson et al, 
2008)) and they may also include stochastic components (Eckermann, 2011). Sub-grid cell 
parameterisations contain parameters that cannot directly be measured in the laboratory 
or in the field (Palmer, 2005), but plausible ranges can be estimated from physical consid- 
erations. Consequently, at least some of these parameters are subsequently 'tuned' so that 
the simulator replicates satisfactorily planetary scale phenomena. For example, they may be 
adjusted to get a satisfactory thermohaline circulation. In this sense, the simulator may be 
viewed as a device that constraints relationships between information formulated at different 
time and spatial levels, namely the grid cell and the planetary scale. In addition to sub-grid 
parameterisations, simulators often involve a number of reasonably ad-hoc assumptions, for 
example about the routing of icebergs in the Atlantic ocean and the amount of freshwater 
subsequently delivered to the ocean surface (Gordon et al, 2000, p. 150). 

There are different purposes to the use of simulators in palaeoclimatology. 

1. One is to evaluate the simulator. The purpose is to show that the simulator reproduces 
palaeoclimate observations satisfactorily (or, at least, more satisfactorily than another 
one), and take this as an element of confidence into future climate predictions with this 
simulator, in the context of the anthropogenic climate change problematic (Braconnot 
et al., 2007, 2012). To this end, experiments are designed carefully to be as realistic 
as possible. 

2. Another purpose is to construct an explanatory framework to a climate phenomenon. 
For example, explaining the desertification of the Sahara 6,000 years ago (Kutzbach 
and Liu, 1997; Claussen et al, 1999), the strong East-Asian monsoon signal during 
marine isotopic stage 13 ( Yin et al, 2008), the timing of the Holocene climate optimum 
(Renssen et al, 2009) or the effect of different oceanic factors to variations in CO2 
(Archer et al, 2000; Chikamoto et al, 2012). Given the nature of simulators, one 
single simulation cannot count as a fully satisfactory explanation of a phenomenon. 
The interest of simulators lies in the possibility of designing series of "intervention" 
experiments (change the forcing, or take the control of certain model components, such 
as vegetation) in order to investigate causal relationships. Even though simulators do 
not replicate reality exactly, these experiments help us to quantify complex causal 
effects, and possibly summarise them as a visual flow diagram (cf. Section 2.1.1). 

3. The third purpose is to generate information that is not immediately accessible from 
palaeoclimate archives. Simulator experiments may be used to provide climate recon- 
structions (Paul and Schafer-Neth, 2005) or constrain unknown quantities, such as the 
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duration of hydrological perturbations associated with Heinrich Events (Roche et al., 
2004). 

The simulator is an imperfect representation of reality and therefore cannot replicate obser- 
vations exactly. It is never 'true' (Oreskes et al., 1994). Consequently, modellers have been 
looking at ways of expressing the distance between the simulation and the reality, either 
qualitatively or quantitatively. The classical procedure is illustrated in Figure 6. The log of 
data produced by the simulation is post-processed and summarised, for example in the form 
of seasonal averages and maps. In parallel, palaeoclimate observations are compiled and 
expressed in a model-friendly (or modeller-friendly!) database, that is, they are expressed in 
terms of climate variables, using statistical or mechanical relationships, and aggregated on 
a grid (an example is the MARGO dataset (MARGO project members, 2009)). Finally, the 
distance between observations and simulation may be discussed either qualitatively (Otto- 
Bliesner et al., 2009), or on the basis of formal metrics. For example, a distance metric 
based on fuzzy logic was proposed in the late 1990's (Guiot et al., 1999). This particular 
metric was designed to avoid penalizing excessively a simulator that would reproduce climate 
change patterns correctly, but with shifts in location compared to reality. The model-data 
comparison procedure is an integral part of the modelling process. More specifically, which 
quantities are being looked at (seasonal averages, regional averages) and which metrics are 
being chosen, reflect the judgements of the modeller about which information generated 
by the simulator is potentially valuable for inference (cf. (Guillemot, 2010) for the science 
historian prospective of on this matter). 

From there research efforts have taken two complementary (and non-exclusive) directions: 

1. Develop process-based models for palaeoclimate observations, so that the simulator 
generates information of which the nature is similar to the one which is observed. The 
process is illustrated on Figure 7. Rather than comparing simulated climate with recon- 
structed climate, one compares simulated isotopes with observed isotopes (Hoffmann 
et al, 1998; Marchal et al, 2000; LeGrande et al, 2006; Roche et al, 2006; Muller 
et al., 2006) or biomes with pollen scores (Harrison and Prentice, 2003; Brewer et al., 
2007). As for climate models, process-based models for palaeoclimate observations are 
available with various levels of sophistication, and they can be intimately embedded 
in the climate simulator (necessary for isotopes simulation) or included as a separate 
post-processing module (e.g.: the BIOME model, Kaplan et al. (2003)). Compared 
to the classical procedure of climate-climate comparison, observation modelling miti- 
gates the loss of information occurring when observations are converted into climatic 
variables. The modeller also has a better control on the assumptions linking climatic 
variables to palaeo-observations. 

2. Develop sampling strategies, founded on ensembles of experiments. For example one 
may sample different parameter combinations of a simulator in order to explore the 
effect of uncertainties on the value of these parameters. The parameter combinations 
giving simulations comparing most favourably to palaeoclimate data may provide more 
trustworthy predictions (Annan et al., 2005; Schneider von Deimling et al., 2006). 
Simulations obtained with these best parameters may also be viewed as a form of 
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Figure 6: Classically, simulations and data-based palaeoclimate reconstructions are developed separately and 
then compared, either based on visual inspection, or using quantitative functions (metrics) measuring the 
distance between the simulation and the palaeoclimate reconstruction. 



climate reconstruction (Paul and Schafer-Neth, 2005). Sampling forcing scenarios may 
also provide constraints on the latter (Roche et ai, 2006). With the ensemble strategy, 
the interest of modellers has moved from the notion of evaluating one experiment to 
evaluating the experiment ensemble, that is, asking whether climate summaries (e.g.: 
zonal averaged temperature) generated by the ensemble is likely to encompass the 
true value (Hargreaves et ai, 2011). Section 3 provides more discussion on how an 
ensemble-based strategy fits a statistical approach. 

2.3. Statistical modelling 

A statistical model is a mathematical construction that accommodates both systematic 
relationships between different variables and random effects (Davison, 2003), and random 
effects are modelled through the notion of probability distributions. 

In classical statistics, one generally distinguishes two objectives: testing a hypothesis, 
and estimating a quantity. In both cases, the inference procedure relies on the definition of 
quantities called 'statistics', and on the specification of a process that could plausibly have 
generated the observations. 

For example, Huybers and Wunsch (2005) ask whether deglaciations occur systematically 
for a same phase of obliquity. To this end they consider a statistic called the "R" statistic, 
and they study the distribution of this quantity that would have been obtained if the process 
that generated the benthic record had been replicated a large number of times, accounting for 
uncertainties on the chronology. They then observe that this hypothesis cannot be rejected 
but another hypothesis, that obliquity phases at terminations are totally random, can be 



15 



simulation of observations 



(uncertain) 
parameters 



Forcing and 
boundary condi- 
tions 



(uncertain) 
parameters 



:limate simulation 



(observation simulation 



summary (e.g. maps) 



paleoobservation 



metric 



Figure 7: As an adaptation of the methodology shown in Figure 6, it is also possible to simulate palaeo- 
observations using appropriate models, for example, vegetation or isotope models, so that the simulation- 
observation comparison is done at the level of the recorded information, rather than at the level of climatic 
variables. 
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rejected. In another application, Braun et al. (2011) use a two-parameter process to test a 
possible solar influence on triggering abrupt events, and they use the 'mean waiting time' 
as the test statistic. As a last example, Haam and Huybers (2010) developed a statistical 
algorithm to test the correlation between two time series with chronological uncertainties. 

So far, most problems of estimation in Quaternary palaeoclimatology were posed as a 
problem of regression, in which it is assumed that an observation is the sum of a systematic 
process and a random component. Namely, observations may be modelled as the sum of a 
climate-controlled component and a random observation error. Depending on the nature of 
the relationship, the regression problem may be solved with more or less complex techniques, 
including linear regression, neural networks and "best analog" approaches (e.g. Kucera et al, 
2005). The problem of interpolating sparse observations on a grid may also be posed as a 
regression problem, and it involves then a model of the spatial dependency of observations. 
Discussion and application on sea-surface temperature estimates of the Last Glacial Maxi- 
mum are available in Schafer-Neth et al. (2005). Finally, specific time-process models have 
been tailored for palaeoclimate time-series analysis. In particular, the 'RAMPFIT' model 
(Mudelsee, 2000) assumes a linear trend between two plateaus. It was used to estimate the 
timing of the Northern Hemisphere glaciation based on foraminifera data (Mudelsee and 
Raymo, 2005). 

More recently, palaeoclimate scientists have become interested in modelling probability 
distributions of uncertain variables, and the Bayesian approach provides a framework to this 
end. The principle is the following. Observations are modelled as the output of a process 
expressed mathematically. The process generally involves a number of parameters. In the 
Bayesian formalism, knowledge on the natural process is being incorporated by different 
channels. The first one is the structure of the equations that represent the process. Examples 
are given later in this review, in which sophisticated vegetation or climate models are used 
to this end; though the modeller may chose a very general model. For example, auto- 
regressive processes are a very broad class of models that express the fact that the system 
has a memory. Physical knowledge also concerns the range of values that parameters may 
plausibly take. This is the second channel. This information is incorporated in the Bayesian 
statistical model under the form of prior probability distributions of parameters, which 
reflect the judgements of the modeller about the probability that the parameters will take 
this or that value. Finally, observations are accounted for in a mathematical procedure that 
consists in updating the prior distributions on parameters by application of the Bayes' law 
of conditional probabilities. This operation is achieved by resorting to a function called the 
likelihood. This is a function of the parameters, which expresses the degree of consistency 
between the model output obtained with these parameters, and observations at hand 7 . The 
resulting probability distribution of the parameters, which combines prior and likelihood, is 
called the posterior probability distribution. 



'Technical definitions of the likelihood and explanation of Bayesian calculus are available in classical 
textbooks, such as Davison (2003), but the broader discussion by (Jaynes, 2003, chap. 3) is also quite 
enlightening. 
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The Bayesian formalism is also equipped with a concept to quantify the degree of self- 
consistency between these different informations. The relevant quantity is called the marginal 
likelihood: Among different models, the one with the highest marginal likelihood is generally 
to be preferred. Though, this is by no means automatic. Assessing a model remains a 
broad and partly subjective procedure that involves testing it against independent data, and 
discussing its structure based on physical considerations. 

Bayesian models were proposed to estimate distributions of sediment age based on un- 
certain radiocarbon dates (Blaauw and Christen, 2005; Haslett and Parnell, 2008), infer 
the posterior distribution of Holocene temperatures based on chironomid taxa abundances 
(Korhola et al, 2002), or estimate ice core chronologies (Lemieux-Dudon et al, 2010). 

The Bayesian approach also provides the flexibility needed to construct 'hierarchical 
models' (Figure 8). The hierarchical model chains different sub-models to form a consistent 
forward description of the process that leads to the observation. For example, it is possible to 
articulate a process to model response of climate to a forcing, with one model that expresses 
the response of a climate record to climate (for example: <5 18 in a marine sediment). The 
statistical resolution of the problem consists in finding the joined probability distribution of 
all variables involved in the model, based on prior evidence and observations. 

Hierarchical modelling was recently introduced in palaeoclimate applications. Haslett 
et al. (2006) use a hierarchical model to reconstruct climate from pollen samples on the 
site of Glendalough in Ireland. Li et al. (2010) apply a hierarchical model to estimate the 
evolution of annual mean, northern hemisphere temperature during the last millennium based 
on three observation time series and three estimates of the climate forcing (volcanoes, solar 
and greenhouse gases), and Tingley and Huybers (2010) combined model for the geographic 
correlation and correlation of climate states in a model, which they then chained with a 
model for the observations. They called this model 'BARCAST', for "Bayesian Algorithm for 
Reconstructing Climate Anomalies in Space and Time" , and they compare it with different 
other approaches, Bayesian or not, in Tingley et al. (2012). 

3. Combining statistical with process-based modelling 

3.1. Principle 

The purpose of process-based modelling is to articulate physical and biogeochemical 
principles, be they deterministic or stochastic. Statistical modelling, on the other hand, is 
concerned with inference about uncertain variables. Here we show how these two approaches 
may be combined. To illustrate the principle, consider the work of Guiot et al. (2000). In this 
study, the BIOME3 process-based vegetation model (Haxeltine and Prentice, 1996) is used 
to generate prior distributions of the state vegetation in modern and past climates. The prior 
is then updated based on observed pollen fossil samples. The approach is thus statistical in 
principle, but physical or biogeochimical principles are used to restrict the range of scenarios 
that may be admitted as plausible. Likewise, Annan et al. (2005) considered a range of 
parameters of the climate simulator MIROC2 to generate a plausible probability distribution 
of the climate of the Last Glacial Maximum. They obtained a posterior distribution of the 



18 



Hierarchical statistical modelling 



climate params 





observation 
params 


forcing 














climate 




observations 



Figure 8: The hierarchical model is a statistical model, which can be represented graphically using a graphical 
network as shown here. The purpose of the network is to represent what each variable depends statistically 
on. In this example, the state of climate depends on parameters (describing climate processes) and on the 
forcing. What is being actually observed (observations) depends on climate and on a number of parameters 
that control the way climate information is being stored in climate's natural archive. Statistical inference 
consists in updating our knowledge of climate and parameters, given observations and knowledge of the 
forcing. 

parameters of that model, which they used to construct a posterior distribution of climate 
sensitivity to anthropogenic perturbations. 

Unfortunately, process-based models are imperfect. Hence, the output of a process-based 
model cannot be taken as a plausible exact representation of reality, even if the model 
parameters are carefully chosen (Goldstein and Rougier (2006) and Rougier (2007)). Conse- 
quently, it was argued that one cannot really be satisfied with sampling model parameters in 
order to generate a prior distribution of plausible climates Rougier (2007). Statisticians and 
climate scientists must therefore co-operate in order to develop a representation of the resid- 
ual distance between the climate simulator and reality. This distance is called discrepancy 
(Goldstein and Rougier, 2006). The discrepancy is a model in itself. It usually involves pa- 
rameters, which can be defined a priori and updated. For example, Schmittner et al. (2011) 
use the combination of a constant bias and a spatially-covariate process to link the results 
of the UVic model to continental and marine palaeodata, although more work is still needed 
to quantify the effect of this discrepancy on current estimates of climate sensitivity. These 
considerations lead us to a proposal programme for inference combining simulators with sta- 
tistical reasoning that follows the schematics shown on Figure 9. At least three challenges 
are being faced in this enterprise: sampling priors, expressing discrepancy, and identifying a 
tractable resolution algorithm. The present review is concerned with sampling priors only, 
because this is the area where progress has been the most tangible over the recent years. 



3.2. Time slice modelling 

To illustrate the challenge consider the problem of reconstructing the climate of the Last 
Glacial Maximum. Suppose the climate-process model is a simulator like HadCM 3 ( Gordon 
et al., 2000), or any similar climate simulator. Experts are being elicited to identify a short- 
list of, say, 10 uncertain input parameters, which are thought to influence the simulated 
climate significantly The problem is to sample this 10-dimension parameter space. If, for 
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Figure 9: Information provided by climate simulators may be involved in a hierarchical model similar to 
Figure 8. Climate is modelled as the combination of climate simulations and discrepancy; and observations 
are modelled as the combination of a process-based observation model (for example: a vegetation model 
such as BIOME) and discrepancy. 
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every parameter, 5 different values are being tried, 5 10 experiments need to be run. That 
makes a bit less that 10 million experiments, each of which takes several weeks of super- 
computer time. This is naturally impossible. 

On this question, the literature opposes two approaches. The first one consists in search- 
ing methodically the region of maximum likelihood (or, maximum posterior) and estimate 
the local curvature of the likelihood (or, posterior) around it. To this end, the modeller may 
resort to the adjoint of the simulator. This is a mathematical or software device that esti- 
mates sensitivity (formally: the gradient) of the model state with respect to its parameters. 
The adjoint method was recently illustrated in the palaeoclimate context (Paul and Losch, 
2012). The reader is referred to textbooks for a more formal definition adjoint approach, 
including Kalnay (2002) and Wunsch (2006). The approach is efficient, but it may fail if 
the posterior has several local maxima. 

The second approach is to sample the 10-dimensional parameter space sparsely but 
smartly. The challenge is then to chose carefully the experiments to be run, such as to extract 
as much information as possible. This is a problem of 'experiment design for computer ex- 
periments', for which a theory is being developed (Santner et ai, 2003), and several articles 
focus explicitly on climate applications (Urban and Fricker, 2010; Challenor, 2011). Once 
experiments have been run, the output of the simulator has to be interpolated throughout 
the entire 10-dimensional space. This reminds us of the spatial interpolation problem (in- 
terpolate sparse observations on a regular mesh). Unsurprisingly, interpolation methods for 
computer experiments use some of the mathematics of spatial modelling and, in particular, 
the Gaussian process (Rasmussen and Williams, 2005). The literature on the subject mainly 
relies on the Bayesian formalism, and the interpolator is called in this context an 'emulator' 
(Kennedy and O'Hagan, 2000; Oakley and O'Hagan, 2002). The modeller's choice may be 
to emulate global quantities, such as the globally averaged temperature (Olson et ai, 2012). 
In practice, though, palaeoclimate applications often require to reconstruct the spatial dis- 
tribution of temperature, precipitation, and other climatic quantities. Several solutions have 
been proposed, including the Principal-Component- Analysis-emulator ( Wilkinson, 2010) and 
the outer-product emulator (Rougier, 2008). Both options were tested in Schmittner et al. 
(2011) (in their supplementary material). 

3.3. Dynamical modelling 

Let us now focus on the dynamical climate evolution over Quaternary time scales. The 
role of the climate process-based model is to provide a physical prior on the physical rela- 
tionships between the different climatic variables. 

Here, one should distinguish two different possible objectives. One objective is to recon- 
struct a coherent picture of climate evolution, using physical constraints. This is a problem 
of state estimation. Another objective is to enquire about climate mechanisms by testing 
hypotheses and estimate parameters. The second task can be translated into a Bayesian sta- 
tistical framework as a problem of parameter estimation. Here we review state and parameter 
estimation separately. 
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3.3.1. State estimation 

In principle state estimation in dynamical systems is a fairly standard problem for which 
a number of solutions exist, based on the Kalman filter (e.g. Evensen and van Leeuwen, 2000) 
or the particle filter (Doucet et ai, 2001). In practice, palaeoclimate applications involve 
numerical systems for which a direct application of theses technologies would be overwhelm- 
ingly expensive. Consequently, somewhat ad hoc, application-dependent algorithms were 
suggested. 

For example, Bintanja et al. (2005) use a thermodynamical ice-sheet model, including an 
explicit representation of oxygen isotope balance within the ice sheet, coupled to a simple 
representation of the deep ocean and the atmosphere. Given observations of only certain 
variables of the model (here: deep ocean 5 18 0), the problem is to reconstruct the state of all 
system variables, and, in particular, deep-ocean temperature and continental ice mass. This 
is a typical problem of state estimation, which was solved in this case by a time-step wise 
iterative procedure. 

Another palaeoclimate application is the spatio-temporal reconstruction — for example, 
reconstructing the geographical distribution of temperature throughout the last millennium — 
using climate simulators. At the palaeoclimate time scale the only tractable approach today 
relies on low resolution general circulation models belonging to the category of climate mod- 
els of intermediate complexity. If the ocean or atmosphere dynamics are chaotic (this is the 
case in 'LOVECLIM', Goosse et al. (2005)) then one needs to assimilate observations se- 
quentially, that is, one after the next as time progresses. Dubinkina et al. (2011) and Goosse 
et al. (2012) have proposed to solve this problem with a low-computational version of the 
particle filter. 

3.3.2. Parameter estimation 

Calibrating the parameters in a model on observations allows one to update prior knowl- 
edge on the relationships between different climate variables. The most tractable solution 
(i.e., the least computational intensive) is to rely on a low-order dynamical system, for ex- 
ample a conceptual model, of which the execution is almost instantaneous on a computer. 
Even though the model is highly idealised, inference on model parameters may provide in- 
formation on causal relationships between the different components of the Earth system. 
For example, if one parameter controls the effect of precession on ocean circulation, it is 
informative to know if this parameter is with high probability different to zero. One model 
proposed in the literature for such applications is the Saltzman and Maasch (1990) model. 
The methodological developments aiming at using this model in a statistical setting have 
paralleled those for snap-shot reconstructions discussed in the previous subsection: Early 
attempts are based on sampling uncertain parameters assuming that the model is perfect 
(Hargreaves and Annan, 2002), and later attempts have included a stochastic component 
to account for model discrepancy (Crucifix and Rougier, 2009). If the stochastic model is 
unstable, i.e., if two climate reconstructions simulated with two different but plausible real- 
isations of the model error term vary widely, then the introduction of stochastic processes 
makes the parameter estimation task a particularly difficult problem, about which litera- 
ture is scattered between the traditions of physics and dynamical system theory (see Voss 
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Figure 10: In the dynamic emulator project, for which Young and Ratto (2011) provide theoretical back- 
ground and applications, the simulator is identified to a reduced, stochastic model, operating in a reduced 
space. For example, instead of the full geographical distribution, the stochastic dynamical system works 
on the space of principal components. For every set of parameters of the simulator, a corresponding set 
of parameters for the stochastic dynamical system can be estimated such as to reproduce results obtained 
with the full simulator as well as possible. The dynamical emulator then consists of the combination of the 
stochastic dynamical system with a statistical model, mapping simulator parameters on stochastic model 
parameters. 

et al. (2004) for a general review, and Kwasniok and Lohmann (2009) for an application on 
Dansgaard-Oeschger events), and statistics (see Ionides et al. (2006); Liu and West (2001) 
and Andrieu et al. (2010) for three different approaches to the problem). 

It may be objected that conceptual models contain too little physical information, and 
that using more comprehensive simulators is desirable. The objection has merit, but how to 
deal with the subsequent intractable computational would unavoidably emerge? Young and 
Ratto (2011) proposed a strategy to address this specific problem, relying on the technique 
of dynamical system identification. This a two-step process (Figure 10). 

The first step consists in finding a stochastic dynamical system that reproduces ad- 
equately the ouptut of a simulator that encapsulates the desired physics. A number of 
practical examples are already available in the literature: 

1. Hooss et al. (2001) considered an ocean carbon cycle model developed to compute the 
oceanic uptake of CO2 called HAMOCC1. HAMOCC1 is a complex simulator, resolv- 
ing the three-dimensional fluid dynamics and a number of oceanic tracers involved in 
the carbon and nutrient cycles. The model outputs three-dimensional fields. Hooss 
et al. (2001) estimated the principal spatial patterns of the simulator output, and cali- 
brated a finite impulse response model predicting the response of the coefficients of the 
principal patterns, in order to reproduce the full output of HAMOCC1. They further 
corrected the model slightly to account for second-order non-linearities related to CO2 
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uptake. 



2. The LLN2D model is a simulator of the growth and decay of the three major northern 
ice sheets, in response to variations in carbon dioxide and astronomical forcing ( Gallee 
et al, 1992). The proposed surrogate of this model (Crucifix, 2011) is a second-order 
(hence, non-linear) equation governing the total amount of continental ice, suitable to 
account for the co-existence of two stable steady-states in this model. 

3&4. Finally, Young and Ratto (2009) considers, among others, two applications: (a) the 
Enting and Lassey (1993) model of the global carbon cycle, which is itself a some- 
what elaborate dynamical system, with 23 coupled nonlinear equations and (b) the 
HadCM3 climate simulator (Gordon et al, 2000), used in response to anthropogenic 
greenhouse gas increases scenarios. For both applications a surrogate model is obtained 
by calibrating a 'Transfer function model' 8 

The second step consists of the Young and Ratto (2011) strategy is to model the relationship 
between the parameters of the simulator and the parameters of the surrogate stochastic 
dynamical system using a Gaussian process. One then obtains a solution to explore the 
response of the simulator to different combinations of parameter and forcing, by actually 
running the surrogate stochastic dynamical system. In this context, Young and Ratto (2011) 
calls the combination of the surrogate dynamical system with the model of the relationship 
between the parameters of the simulator and those of the surrogate system, a dynamical 
emulator. So far no palaeoclimate application of this technique is available. 

3.4- Emulator by-products 

The emulator was introduced and defined as a statistical surrogate of a simulator used 
to sample efficiently the space of uncertain parameters. More generally, the whole process 
of designing and calibrating an emulator is a procedure aiming at extracting information 
efficiently from a complex numerical system. This is a computing cost economiser, which 
comes with at least two possible applications in addition to the problems of parameter and 
state estimation discussed above: 

1. Once the emulator is properly calibrated and tested, there is almost no limit to the 
number of sensitivity experiments to model inputs that may be run, because they 
are run with the emulator rather than with the actual simulator. Consequently, the 
sensitivity of the simulator can be explored very finely and span the whole input space, 
which includes boundary conditions, forcing and parameters. The procedure provides 
us with a way to characterise the non-linear structure of the model response and, in 
particular, identify thresholds and plateaus. Oakley and O'Hagan (2004) provide some 
theoretical background to this end. 



8 Transfer-function models are classical in dynamical system identification theory (Ljung, 1999). The 
assumption is that moving-averaged outputs and moving averaged- inputs are linearly correlated to each 
other. 



24 



2. Given the computational resources required to run general circulation models in palaeo- 
climate applications, modellers generally chose to degrade the spatial resolution and 
simplify the physics (for example: assume that cloud cover is constant, or simplify 
the relationship between large scale atmospheric dynamics and precipitation) com- 
pared to state of the art simulators. This is the rationale for the use of Earth System 
Models of Intermediate Complexity (Claussen et ai, 2002). However, these physical 
simplifications have deleterious effects on the representation and the response struc- 
ture of certain processes thought to play a role at the global scale, including tropical 
monsoons, El-Nino, Southern ocean ventilation and North Atlantic convection. The 
emulator provides an alternative to this reduced-physics strategy, by offering the pos- 
sibility of running numerous experiments with an emulator calibrated on a simulator, 
rather than with the simulator itself. 

This being said, developing an emulator is by no means straightforward, be it in the static or 
in the dynamic framework. The choice of emulator structure, inputs, and regressors is a fairly 
heuristic procedure that requires expertise on the simulator and the processes that it is meant 
to represent. In particular, one cannot overstate the importance of extensive validation tests 
of the emulator: the user must be confident that the emulator is consistent with the simulator 
dynamics and avoid unsuitable interpretations. The literature on emulator validation and 
design is developing (Bastos and O'Hagan, 2009) (cf. also the www.mucm.ac.uk website) and 
should be carefully examined before any application. The natural recommendation, echoed 
by Young and Ratto (2011) (and references therein) is to adopt a framework that is as simple 
as possible, chosen such that the parameters of the surrogate model can be interpreted in 
terms of system dynamics and system sensitivity. 

4. Conclusion 

The present article reviewed three modelling approaches adequate to study palaeoclimate 
phenomena: 

1. Conceptual models provide a representation of a well-identified phenomenon like ice 
ages, Dansgaard-Oeschger events or Heinrich events. These models are specific to the 
studied phenomenon, and their development obeys a principle of parsimony: explain- 
ing as much as possible with as little as possible. In spite of this principle of parsimony, 
there may be a risk of over-interpretation. Two conceptual models relying on different 
mechanistic interpretations may have similar mathematical descriptions, and thus fit 
the same dataset equally well (Tziperman et ai, 2006; Cane et ai, 2006; Crucifix, 
2012). Consequently, a model fit does not necessarily imply that the underlying causal 
interpretation is correct. As a general rule, the model will be more convincing if it 
accommodates different forms of information, that is: different records (5 18 0, CO2, 
SST,. . . ), different aspects of the record (trends, noise, auto-correlation,. . . ), and cap- 
tures transitions between different regimes (for example: the Middle Pleistocene Revo- 
lution) with few free parameters. The theory of dynamical system tells us that chaotic 
and other counter-intuitive mathematical phenomena may arise in deceptively simple 
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systems. For this reason, dynamical system identification, which we briefly reviewed at 
the end of the article, is useful as a tool for analysing the dynamics of simulators used 
to generate time series. If the behaviour of a complex simulator can be reproduced 
with a compact, stochastic dynamical system, it can be more easily analysed from the 
prospective of dynamical system theory. 

2. Simulators provide a spatially-resolved representation of numerous physical and bio- 
geochemical processes involved in the climate system. Simulators are more versatile 
in nature than conceptual models. Depending on the experimental design, the same 
climate simulator may be used to study different climatic phenomena at different time 
and/or spatial scales. Causal effects may be explored by means of intervention exper- 
iments, in which one of the components of the system is being taken control of (e.g.: 
fixed or 'imposed' ice sheets, vegetation, CO2 etc.) However, it would be a mistake to 
regard simulators as ab initio models of the climate system. To the extent simulators 
and their sub-grid scale parameterisations are tuned to adequately reproduce planetary 
scale phenomena (globally-average temperature, the Atlantic meridional overturning 
cell, monsoon dynamics), simulators are partly empirical. Scientists working with cli- 
mate simulators are aware of this and for this reason they develop stringent tests to 
assess the performance of these models, both for the representation of the present-day 
climate and for the past climates. 

3. Statistical modelling accommodates systematic relationships — but not necessarily causal— 
between different variables with random effects. Statistical modelling may in particular 
be involved in space-time climate reconstruction, hypothesis testing and parameter es- 
timation, but it should never be overlooked that these mathematical operations are 
constrained by a number of hypotheses expressed in the statistical model. Common 
statistical assumptions concern the autocorrelation of data in time and space, and 
between the climate state and actual observations. Hierarchical modelling is increas- 
ingly used because it provides a more explicit representation of the forward mechanical 
process linking the causes of the phenomenon to actual observations. 

One of the purposes of the review is to show that statistics, and Bayesian statistics in 
particular, have a role to play to articulate these frameworks. In particular, they allow us to 
include information on a priori physically plausible model parameter ranges in the inference 
process. 

Unfortunately there is no free lunch. A brute-force sampling approach may involve 
millions or even billions of experiments, which cannot be envisaged even with the simplest 
systems. Consequently, resolution algorithms tend to be mathematically highly involved, 
and resort to intermediate modelling steps such as the emulator or the adjoint, experiment 
planning, along with new and not yet solidly established concepts such as the discrepancy. 
These modelling steps and associated resolution algorithms need to be evaluated carefully 
and this procedure is tedious and technical. What has been gained on the one hand (a 
more transparent expression of the conditions of inference) has to be paid by the technical 
character of the inference process, which may appear obscure and off-putting to a significant 
fraction of the scientific audience. 
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The conclusion follows naturally. The established paradigms of modelling that were re- 
viewed here (conceptual, simulator and statistical modelling) still have an important role to 
play because they allow scientists to build on years of achievement and experience. Concep- 
tual modelling, in particular, is a flexible approach to organise observations involved in com- 
plex climatic processes. On the other hand, the power of cross-paradigm approaches, allying 
statistics with simulators, or statistics with conceptual models, should not be overlooked. 
They offer avenues to manage uncertainties and rationalise model design and selection given 
different and variate sources of information. 
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